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Abstract 

A time-domain numerical modeling of transversely isotropic Biot poroelastic waves is proposed 
in two dimensions. The viscous dissipation occurring in the pores is described using the dynamic 
permeability model developed by Johnson-Koplik-Dashen (JKD). Some of the coefficients in the 
Biot-JKD model are proportional to the square root of the frequency. In the time-domain, these 
coefficients introduce shifted fractional derivatives of order 1 /2, involving a convolution prod¬ 
uct. Based on a diffusive representation, the convolution kernel is replaced by a finite number 
of memory variables that satisfy local-in-time ordinary differential equations, resulting in the 
Biot-DA (diffusive approximation) model. The properties of both the Biot-JKD and the Biot-DA 
model are analyzed: hyperbolicity, decrease of energy, dispersion. To determine the coefficients 
of the diffusive approximation, two approaches are analyzed: Gaussian quadratures and opti¬ 
mization methods in the frequency range of interest. The nonlinear optimization is shown to be 
the better way of determination. A splitting strategy is then applied to approximate numerically 
the Biot-DA equations. The propagative part is discretized using a fourth-order ADER scheme 
on a Cartesian grid, whereas the diffusive part is solved exactly. An immersed interface method 
is implemented to take into account heterogeneous media on a Cartesian grid and to discretize 
the jump conditions at interfaces. Numerical experiments are presented. Comparisons with an¬ 
alytical solutions show the efficiency and the accuracy of the approach, and some numerical 
experiments are performed to investigate wave phenomena in complex media, such as multiple 
scattering across a set of random scatterers. 
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1. Introduction 


A porous medium consists of a solid matrix saturated with a fluid that circulates freely 
through the pores Gill! . Such media are involved in many applications, modeling for instance 
natural rocks, engineering composites [ 0 ] and biological materials [5]. The most widely-used 
model describing the propagation of mechanical waves in porous media has been proposed by 
Biot in 1956 |' 1, [qD . It includes two classical waves (one ’’fast” compressional wave and one 
shear wave), in addition to a second ’’slow” compressional wave, which is highly dependent on 


the saturating fluid. This slow wave was observed experimentally in 1980 00, thus confirming 
the validity of Biot’s theory. 

Two frequency regimes have to be distinguished when dealing with poroelastic waves. In 
the low-frequency range (LF), the flow inside the pores is of Poiseuille type (3]. The viscous 
effects are then proportional to the relative velocity of the motion between the fluid and the 
solid components. In the high-frequency range (HF), modeling the dissipation is a more delicate 
task. Biot first presented an expression for particular pore geometries fl]. In 1987, Johnson- 
Koplik-Dashen (JKD) published a general expression for the HF dissipation in the case of random 
pores HI, where the viscous efforts depend on the square root of the frequency. No particular 
difficulties are raised by the HF regime if the solution is computed in the space-frequency domain 
mu. On the contrary, the computation of HF waves in the space-time domain is much more 
challenging. Time fractional derivatives are then introduced, involving convolution products 
0. The past of the solution must be stored, which dramatically increases the computational 
cost of the simulations. 

The present work is proposes an efficient numerical model to simulate the transient poroe¬ 
lastic waves in the full frequency range of Biot’s model. In the high-frequency range, only two 
numerical approaches have been proposed in the literature to integrate the Biot-JKD equations 
directly in the time-domain. The first approach consists in a straightforward discretization of 
the fractional derivatives defined by a convolution product in time Gi. In the example given 
in 11211 . the solution is stored over 20 time step s. The second approach is based on the diffusive 
representation of the fractional derivative lfl3n . The convolution product is replaced by a con¬ 
tinuum of memory variables satisfying local differential equations Gi. This continuum is then 


discretized using Gaussian quadrature formulae | 


17H . resulting in the Biot-DA (diffusive 


approximation) model. In the example proposed in [ 13fl, 25 memory variables are used, which is 
equivalent, in terms of memory requirement, to storing 25 time steps. The idea of using memory 
variables to avoid convolution products is close to the strategy commonly used in viscoelasticity 

Q. 


The concern of realism leads us also to tackle with anisotropic porous media. Transverse 
isotropy is commonly used in practice. It is often induced by Backus averaging, which re¬ 
places isotropic layers much thinner than the wavelength by a homogeneous isotropic transverse 
medium Cl. To our knowledge, the earliest numerical work combining low-frequency Biot’s 
model and transverse isotropy is based on an operator splitting in conjonction with a Fourier 
pseudospectral method [ 20]. Recently, a Cartesian-grid finite volume method has been developed 
[21]. One of the first work combining anistropic media and high-frequency range is proposed in 
1 2211 . However, the diffusive approximation proposed in the latter article has three limitations. 
Firstly, the quadrature formulae make the convergence towards the original fractional operator 
very slow. Secondly, in the case of low frequencies, the Biot-DA model does not converge to¬ 
wards the Biot-LF model. Lastly, the number of memory variables required for a given accuracy 
is not specified. 
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The present work extends and improves our previous contributions about the modeling of 
poroelastic waves. In ll23ll . we addressed ID equations in the low-frequency range, introducing 
a splitting of the PDE. 2D generalizations for isotropic media required to implement space-time 
mesh refinement [24, 25]. Diffusive approximation of the fractional derivatives in the high- 
frequency range were introduced in ]|26h and generalized in 2D in [[ 27 I 1 . Compared with [ 27], the 
originality of the present paper is threefold: 


1 . 


incorporation of anisotropy. The numerical scheme and the discretization of the interfaces 
need to be largely modified accordingly; 


2. new procedure to determine the coefficients of the diffusive approximation. In 11261 12711. 


we used a classical least-squares optimization. It is much more accurate than the Gauss- 
Laguerre technique proposed in Ill3h . But in counterpart, some coefficients are negative, 
which prevents to conclude about well-posedness of the diffusive model. Here, we fix this 
problem by using optimization with constraint of positivity, based on Shor’s algorithm. 
Moreover, the accuracy of this new method is largely improved compared with the linear 
optimization; 

3. theoretical analysis. A new result about the eigenvalues of the diffusion matrix is intro¬ 
duced and the energy analysis is extended to anisotropy. 


This article is organized as follows. The original Biot-JKD model is outlined in § [2] and the 
diffusive representation of fractional derivatives is described. The energy decrease is proven, 
and a dispersion analysis is done. In § 0 an approximation of the diffusive model is presented, 
leading to the Biot-DA system. The properties of this system are also analyzed: energy, hyper- 
bolicity and dispersion. Determination of the quadrature coefficients involved in the Biot-DA 
model are investigated in § 13.41 Gaussian quadrature formulae and optimization methods are 
successively proposed and compared, the latter being finally preferred. The numerical model¬ 
ing of the Biot-DA system is addressed in § 01 where the equations of evolution are split into 
two parts: the propagative part is discretized using a fourth-order finite-difference scheme, and 
the diffusive part is solved exactly. An immersed interface method is implemented to account 
for the jump conditions and for the geometry of the interfaces on a Cartesian grid when dealing 
with heterogeneous media. Numerous numerical experiments are presented in §[5] validating the 
method developed in this paper. In §0 a conclusion is drawn and some futures lines of research 
are suggested. 


2. Physical modeling 


2.1. Biot model 

We consider a transversely isotropic porous medium, consisting of a solid matrix saturated 
with a fluid that circulates freely through the pores d&i . The subscripts 1 , 3 represent the x, z 
axes, where z is the symmetry axis (figure[T]). The perturbations propagate with a wavelength A. 

The Biot model involves 15 positive physical parameters: the density p/, the dynamic vis¬ 
cosity rj and the bulk modulus Kf of the fluid, the density p s and the bulk modulus K s of the 
grains, the porosity 0 ^ 0 ^ 1 , the tortuosities Tj > 1 , T 3 ^ 1 , the absolute permeabilities at 
null frequency k\ , K 3 , and the symmetric definite positive drained elastic matrix C 
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Figure 1: medium under study, (a): the physical properties are symmetric about the axis z that is normal to the plane 
(x,y) of isotropy, (b): interface T separating two poroelastic media Ho an d Hi. The normal and tangential vectors at a 
point P along T are denoted by n and t, respectively. 


The linear Biot model is valid if the following hypotheses are satisfied [28]: 


*H\ : the wavelength A is large in comparison with the characteristic radius of the pores r; 

7^2 : the amplitudes of the waves in the solid and in the fluid are small; 

7^3 : the single fluid phase is continuous; 

^4 : the solid matrix is purely elastic; 

7^5 : the thermo-mechanical effects are neglected, which is justified when the saturating fluid 
is a liquid. 


In the validity domain of homogenization theory (TTi), two frequency ranges have to be distin¬ 
guished. The frontier between the low-frequency (LF) range and the high-frequency (HF) range 
is reached when the viscous efforts are similar to the inertial effects. The frequency transitions 
are given by fl] 


/a = =^- = £4 i = 1,3. 

2 nTiKipf 2 n 


( 2 ) 


Denoting u s and Uf the solid and fluid displacements, the unknowns in a velocity-stress for¬ 
mulation are the solid velocity v s = ^ 7 , the filtration velocity w = = f- f (p (uf - u s ), the 

elastic symmetric stress tensor a and the acoustic pressure p. Under the hypothesis of small 
perturbations OH 2 ), the symmetric strain tensor s is 


e = 1 (Vm s + Vh/). (3) 

Using the Voigt notation, the stress tensor and the strain tensor are arranged into vectors cr and s 
a = (cm , 0 - 33 , 0 -i 3 ) r , e = (£ 11 , £ 33 , 2 si3) T . (4) 
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Setting 


£, = -V.'W, C u =C + mPP T , 


(5a) 


a tn q a \T o i C 11 + c 12 + Cl3 , 2 C 13 +C 33 

P = (Pi, Pi, Pi) 1 , p i = l-— -, Pi = 1-—-, (5b) 


3 K, 


3 K, 


K = K s (l+4>(K s /K f -l)), m 


Ki 


K — (2 c\\ + c 33 + 2 ci 2 + 4cn)/9 


(5c) 


where C u is the undrained elastic matrix and 4 the rate of fluid flow, the poroelastic linear con¬ 
stitutive laws are 0 

cr = C“ s - in p 4, p = m (f ~p T s). ( 6 ) 

Using (l5ab and (l5bk we obtain equivalently 

cr = Cs-J3p , p = m(^-J3 r s}. (7) 

The symmetry of a implies compatibility conditions between spatial derivatives of the stresses 
and the pressure, leading to the Beltrami-Michell equation [291120] 


(9Vi 3 ^ d 2 cr\\ ^ d 2 cr 33 ^ d 2 p ^ <9 2 <xn 


n n = 

dxdz 

©0 = -^55 ■ 


d x 2 

C\3 


+ 0i 


n 9 -1" 02 ~Z 9 + ©3 n 9 

or ox 1 oz 


^ d 2 cr 33 d 2 p 

+ 0 o— + 0 <rf’ 


C 11 C 33 -c 


0 i 


13 


£ll 

^13 


00, 02 = /?l 00 + /?3 ©1, 


( 8 ) 


C 33 


03 =- 0 o, 04 = /?3 0 o + /?1 ©3- 

^13 

If the medium is isotropic and in the elastic limit case (J3\ = - 0), we recover the usual 

equation of Barre de Saint-Venant. 

Introducing the densities 


7“ 

P = 0P/ + (l-0)Ps, p w i — ~T Pf 9 i= 1,3, 


(9) 


the conservation of momentum yields 


dv s dw _ 
p - 7 ^- + pf - 77 - = V . <T, 


dt 


dt 


Pf 


dv s dw In \ 

— + diag (p^) — + diag - F;(f) * w = -V p, 


dt 


dt 


\Ki 


( 10 a) 

( 10 b) 


where diag(d ; ) denotes the 2 x 2 diagonal matrix ° ), * denotes the time convolution prod¬ 
uct and F/(t) are viscous operators. In LF, the flow in the pores is of Poiseuille type, and the 
dissipation efforts in (II Obi) are given by 


Fi(t) = F FF (t) = 5{t) <=> Ff (0 * wfa, z , t) = wfa, z , 0, f = 1,3, 




(ID 


where c> is the Dirac distribution, which amounts to the Darcy’s law. 
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2.2. High frequency dissipation: the JKD model 


In HF, a Prandtl boundary layer occurs at the surface of the pores, where the effects of vis¬ 
cosity are significant. Its width is inversely proportional to the square root of the frequency. Biot 
(1956) presented an expression of the dissipation process for particular pore geometries |@]. A 
general expression for the viscous operator for random networks of pores with constant radii has 
been proposed by Johnson, Koplik and Dashen (1987) [ 8 ]. This function is the most-simple one 
fitting the LF and HF limits and leading to a causal model. The only additional parameters are 
the viscous characteristic length A* . We take [12] 


_ 4 7~i Ki n = (Oci = i)(f> 2 A 2 

' </> Af’ ' Pi 47? k 2 


( 12 ) 


where P, is the Pride number. The Pride number describes the geometry of the pores: P ; =1/2 
corresponds to a set of non-intersecting canted tubes, whereas P, = 1/3 describes a set of canted 
slabs of fluids [3]. Based on the Fourier transform in time, Ffoj) = T (Fft)) = L Fft)e~ jojt dt , 

i—I JK. 

the viscous operators given by the JKD model are [ 8 ] 


( 47? K}p f \ m ( u>\ l/2 1 

i+./^ * '7 = l+jp,— =-=(a 

v /;A^0 2 ) \ CL> ci ] 


F] {(F) = \l + j(L> 

Therefore, the terms F ft) * wfx , z, t) involved in (TfObt are 


+ jcj) 1/2 . (13) 


Fi KD (t) * Wi (x, z, t) = T~ x |-2= (Q, + ja» ll2 Wi{x, z, w)J, 


= —j= (D + Q.i) l/2 Wi(x, z, t). 


(14) 


In the last relation of <m, (D + D /) 1/2 is an operator. D 1/2 is a fractional derivative in time of 
order 1/2, generalizing the usual derivative characterized by ^7 = ( jajWi ). The notation 

C D + Q ;) 1/2 accounts for the shift Q, in (IT4l) . 


2.3. The Biot-JKD equations of evolution 


The system (ITOl) is rearranged by separating ( ^f and 
Taking 


- — —L_ i = i 

^ Xi vzr 


and using the definitions of s and f. 
3, (15) 
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one obtains the following system of evolution equations 




3v s i 

Pwl 

3? 

Ti 

3 v s3 

Pw3 

d£ 

X3 

d Wi 

+ eii 

d£ 

Xi ' 


/dan 

da u 

\_P£^£ 

\ dx 

dz 

I Xi dx 

/dcri 3 

da 33 

\_P£?_P 

\ dx 

dz 

I X3 dz 

dan ! 3 < t 13 \ 

+ P_dp__ 

dx 1 dz j 

Xi dx 


= — n (D + ^i ) 1/2 wi + G Vl , 
p 

--^y 3 (D + Q 3 ) 1/2 w 3 +G Vs3 , 
—yi (D + f^i) 1 ^ 2 Wi + G Wl , 


dw 3 pf ( do-13 do-33 \ P dp 

dt + X3 \ dx + dz ) + X 3 dz 


-73 iP + ^ 3) 1/2 ^3 + G W3 , 


d o~n 
d£ 


- c 


11 


dvji 

dx 


- c 


13 


dv,3 

dz 


-mpi 


dwi dw 3 \ 

3x + 3z ]~ G<r “’ 


dcr 13 „ (dv s3 1 3v*i\ _ ^ 

3? 55 (3i 3z ] _tjff13 


do -33 _ „ 3v,i _ u dv , 3 

3? dx ^ 33 Q% 


-mfc 


dw x dw 3 \ 
dx + dz )~ G<t33 ’ 


Ti +m ^' 


3y^i 

dx 


dv s 3 3wi 3 >V 3 

+ “a— + "a- h “3— 

3z 3x dz 


= G„ 


(16a) 

(16b) 

(16c) 

(16d) 

(16e) 

(16f) 

(16g) 

(16h) 


The source terms G Vjl , G Vj3 , G Wl , G w ,, (7 ir|l , G, T|) , G, T) , and G p have been introduced to model 
the forcing. 


2.4. The diffusive representation 

The shifted fractional derivatives in (ITdli can be written l3lll 


(D + Q,) 1/2 w,(x, z, t) 



e -Q.i(t-T) 

ffn (t - r) 


(^y-U z, t) + Q, w,(x, z, r) 


dr, 


i = 1,3. 


(17) 


The operators (Z) + D;) 1/2 are not local in time and involve the entire time history of w. Based 
on Euler’s Gamma function, the diffusive representation of the totally monotone function -j= is 



1 


A JjTt 



-X e~ et dO. 

Xe 


(18) 


Substituting (fl8l) into ([171) gives 


i r°° 1 

(D + £li) l/2 Wi(x , z,t) = - I — z, 6 ,0 dO , (19) 

^ Jo Vd 


where the memory variables are defined as 


if/i(x, z, 0, t) 



e -(e+QiXt-r) 


I dwi 

I Z> T ) + Q W i ( X ’ ^ T ) 
7 


dr. 


( 20 ) 


























For the sake of clarity, the dependence on Eli and w* are omitted in 0 ^. From (l20k it follows that 
the two memory variables if/i satisfy the ordinary differential equation 

( d il/i d Wi 

-qJ = -(0 + EIJ if/i + -qJ- + El t Wi, (21a) 

0 ) = 0 . ( 21 b) 

The diffusive representation therefore transforms a non-local problem (IT 71 ) into a continuum of 
local problems (fT9l) . It should be emphasized at this point that no approximation have been made 
up to now. The computational advantages of the diffusive representation will be seen in §[3]and 
[5J where the discretization of (IT9l) and (12 1 al) will yield a numerically tractable formulation. 


2.5. Energy ofBiot-JKD 

Now, we express the energy of the Biot-JKD model (IT 6 b . 

Proposition 1 (Decrease of the energy). Let us consider the Biot-JKD model f fTbl) without forc¬ 
ing, and let us denote 

E = E l +E 2 + E 3 , (22) 

with 

E\ = - I (pv s T v s + 2pfV s T w +w T diag (p w J w) dxdz , 

2 Jr 2 v ' 

E 2 = ^ f |(<T + p/?) r C- 1 (<r + p/?) + - p 2 ] (iv Jz, (23) 

2 J R 2 \ m / 

£3 — — f - f (w-^) r diag(- -) (w - i//)d6 dx dz. 

2 Jr2 7r J o \ Ki 6 (6 + 2 Q/) / 

E is an energy which satisfies 


dE 

dt 


fir 

Jr 2 n Jo 


ij/ T diag 


6 + flj 


<i s/El; 6 {6 + 2 EIJ 


* 


+w T diag 


Q i 


si Eli 6 {6 + 2 EIJ 


w > d6dxdz ^ 0 . 


(24) 


Proposition \l \ is proven in | Appendix A It calls for the following comments: 


• the Biot-JKD model is stable; 


• when the viscosity of the saturating fluid is neglected (77 = 0 ), the energy of the system is 
conserved; 

• the terms E\ and E 2 in (l23l) have a clear physical significance: E\ is the kinetic energy, 
and E 2 is the strain energy; 

• the energy analysis is valid for continuously variable parameters. 
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2.6. Dispersion analysis 

In this section, we derive the dispersion relation of the waves which propagate in a poroelastic 
medium. This relation describes the frequency dependence of phase velocities and attenuations 
of waves. For this purpose, we search for a general plane wave solution of (O 

I V = (Vi, v 3 , Wi, w 3 ) t = Vo e j(to, - k r \ 

\ T = (cr u , u 13 , cr 33 , -p) T = T 0 e j(cJt ~ k r) , 

where k = &(cos(<£>), sin(<£>)) r is the wavevector, k is the wavenumber, Vo and To are the polar¬ 
izations, r = (x, z) T is the position, oj = 2 n f is the angular frequency and / is the frequency. 
By substituting equation (1231) into equations (116el) - (116hl) , we obtain the 4 x 4 linear system: 


r> u 

C 11 

C \ 3 S(p 

f3 \ mcy 

J3\ m s Ci 

c 55 ‘V 

n 

c 55 L( P 

0 

0 

r u r 
^ 13 ^V 

c 33 

fomcy 

m s { 

Pimcy 

/S 3 m s v 

mcy 

msy 


c 


(26) 


where c v = cos(p) and = sin(</>)- Then, substituting d25b into (1 1 6at >- d 1 6db gives another 4x4 
linear system: 


—k 


' s v 0 O' 

0 c f s<p 0 

0 0 0 

,0 0 0 s v , 


£ 


T = oj 


' p 0 

0 p 

pf 0 

0 pf 


pf 

0 

Y( KD ((o) 

ju 

0 

r 


0 

Pf 

o 

Yl KD (aj) 

jco ) 


V, 


(27) 


where Y J ™ and Y^ KD are the viscodynamic operators [ 32]: 

Y? kd = j Cjp wi + 1 Ff KD (co), i= 1,3. (28) 

Ki 

Since the matrix T is invertible, the equations (f26b and (1371) lead to the eigenproblem 

T l £CV = (|) 2 V. (29) 

The equation ([23b is solved numerically, leading to two quasi-compressional waves denoted qPf 
(fast) and qP s (slow), and to one quasi-shear wave denoted qS 12011 . The wavenumbers thus 
obtained depend on the frequency and on the angle p. One of the eigenvalues is zero with 
multiplicity two, and the other non-zero eigenvalues correspond to the wave modes ±k p f(a>, p), 
±k ps (cj, <p) and ±k s (aj, p). Therefore three waves propagates symmetrically along the directions 
cos(^) v + sin(<£>) z and - cos(^) v - sin(^) z. 
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The wavenumbers give the phase velocities c p f(a >, cp) = oj/%Q(k p f), c ps (oj , p) - col%Q{k ps ), 
and c s {(jj,(p) = (ol%z(k s \ with 0 < c ps < c p f and 0 < c s . The attenuations a p f(co,(p) = 
-3m (k p f), a ps (co,(p) = -3m (k ps ) and a s (oj, p) = -3m (k s ) are also deduced. Both the phase 
velocities and the attenuations of Biot-LF and Biot-JKD are strictly increasing functions of the 
frequency. The high-frequency limits (oj —> oo in d29b ) of phase velocities c™ s (cp) and 

c™(cp) are obtained by diagonalizing the left-hand side of (O. 

Various authors have illustrated the effect of the JKD correction on the phase velocity and on 
the attenuation Q. In figure [2l the physical parameters are those of medium Qo (cf table [B, 
where the frequencies of transition are f c \ - 25.5 kHz, f c 3 = 85 kHz. The dispersion curves are 
shown in terms of the frequency at ip = 0 rad. The high-frequency limit of the phase velocities 
of the quasi-compressional waves are c^(0) = 5244 m/s and c^(0) = 975 m/s, which justifies 
the denomination ’’fast” and ’’slow”. Figure[2| calls for the following comments 


• when / < f ci , the Biot-JKD and Biot-LF dispersion curves are very similar as might be 
expected, since F J i KD { 0) = F LjF (0) = 1; 


• the frequency evolution of the phase velocity and of the attenuation is radically different 
for the three waves, whatever the chosen model (LF or JKD): the effect of viscous losses 
is negligible on the fast wave, small on the shear wave, whereas it is very important on the 
slow wave; 

• when / <$c f ci , the slow compressional wave is almost static H0|. When / > f ci , the 
slow wave propagates but is greatly attenuated. 


Taking 


17, = 


' 1 

0 

0 

0 ' 


' 0 

0 

1 

0 ^ 

0 

0 

1 

0 

, u 3 = 

0 

1 

0 

0 

0 

0 

0 

1 


0 

0 

0 

0 

5 0 

0 

0 

0 > 


V 0 

0 

0 

1 , 


the energy velocity vector V e is [ 2j/ 351: 

(P) 


V e = 


<£> 
c e y 


(E s + Ef) 

(P) = -I tRe((?l(C/i.r) r + el(U 3 .T) T )y), 
(co/k) 2 ^ 


{E) 


1 


Ke 1 + 


\co/k\ 2 


v r rv|, 


(30) 


(31) 


where V is the complex conjugate of V, (P) is the Umov-Poynting vector, ( Ek) and (E s ) are the 
average kinetic and strain energy densities, and (E) is the mean energy density. The theoreti¬ 
cal wavefronts are the locus of the end of energy velocity vector V e multiplied by the time of 
propagation. We will use this property in §0to validate the simulations. 


3. The Biot-DA (diffusive approximation) model 

The aim of this section is to approximate the Biot-JKD model, using a numerically tractable 
approach. 
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phase velocity of the Pf wave 


attenuation of the Pf wave 




phase velocity of the S wave 



attenuation of the S wave 



phase velocity of the P s wave 


attenuation of the P s wave 




Figure 2: dispersion curves in terms of the frequency. Comparison between Biot-LF and Biot-JKD models at (p = 0 
rad. The vertical dotted line denotes the critical frequency separating low-frequency and high frequency regimes. The 
horizontal dotted lines in the left row denote the maximal phase velocity at infinite frequency. 
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3.1. Diffusive approximation 

The diffusive representation of fractional derivatives (IT9l) is approximated by using a quadra¬ 
ture formula on N points, with weights a l £ and abcissae 6 l { (i = 1,3): 


(.D + £li) l/ 2 Wi(x , z, t ) 



1 

Te 


N 

i// l (x, z, 0, t) d6 - ^ a\ i//\x, z, 0\, t), 

i=\ 


= ^a l f i// l fix,z, t). 

t=\ 


(32) 


From (121aF the 2 N memory variables if/ l { satisfy the ordinary differential equations 

dffn dwi 

-7— = -($ + Q) + -7— + Q Wj, 

< at 1 L at 

i// l fix, z, 0) = 0. 


(33) 


3.2. 77ze Biot-DA first-order system 

The fractional derivatives involved in the Biot-JKD system (IT 6 l) are replaced by their diffusive 
approximation (l32k with evolution equations (l33t . After some algebraic operations, the Biot-DA 
system is written as a first-order system in time and in space, used in the numerical simulations 
Of §0(7= 


dv sX 


_ Pwi_ I doji doj 3 \ Pf d_p _ Pf 

dt x\ \ dx dz '' 

dv S 3 _ Pw3 (da 13 do-33 
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Taking the vector of unknowns 


U = 0*1 , V s3 , Wi , W 3 , <x n , cr 13 , cr 33 , p, \f/\ , , ■ ■ ■ , p N , p N ) T , (35) 

and the forcing 

G — (G Vsl , G Vs3 , G W] , G W3 , G(j n , Go- 13 , Gq - 33 , Gp , G W] , G W3 , G W] , G W3 ) , (36) 

the system (l34l) is written in the form: 

?f-+A^+B?f = -SU + G, (37) 

ot ox oz 

where A and B are the (2 N + 8) x (2 N + 8) propagation matrices and S' is the diffusive matrix 
(given in |Appendix B ). The number of unknowns increases linearly with the number of memory 
variables. Only the matrix S depends on the coefficients of the diffusive approximation. 


3.3. Properties 

Some properties are stated to characterize the first-order differential system (l34l ) . First, one 
notes that the only difference between the Biot-LF model, the Biot-JKD model and the Biot-DA 
model occurs in the viscous operators 


Fi(co) 


' F\ f {cj) = 1 

Ff KD (oj) = -k (Q, + jco) l/2 

yO/ 


Biot-LF, 

Biot-JKD, 


(38) 


F? A (oj) = 


Qi + j oj 


z 

€=1 


O l £ + + j oj 


Biot-DA. 


The dispersion analysis of the Biot-DA model is obtained by replacing the viscous operators 
Ff KD (oj) by F® a (oj) in (l28l) . One of the eigenvalues of T _1 X C (l29l) is still zero with multiplicity 
two, and the other non-zero eigenvalues correspond to the wave modes ±k p f(cj, cp), ±k ps (a), ip) 
and ±k s (oj, ip). Consequently, the diffusive approximation does not introduce spurious wave. 


Proposition 2. The eigenvalues of the matrix M = cos(<£>) A + sin(^) B are 

sp(M) = {0 , ±c%(<p), ±c;», ±c»), (39) 

with 0 being of multiplicity 2 N + 2. 

The non-zero eigenvalues do not depend on the viscous operators Ffoj). Consequently, the high- 
frequency limits of the phase velocities c^(tp), c™ s {ip) and defined in § 12.61 are the same 

for both Biot-LF, Biot-JKD and Biot-DA models. An argumentation similar to Il2lll shows that 
the matrix M is diagonalizable for all ip in [0, 2n[, with real eigenvalues. The three models are 
therefore hyperbolic. 

Proposition 3 (Decrease of the energy). An energy analysis of ft34\) is performed. Let us con¬ 
sider the Biot-DA model f f7?l) without forcing, and let us denote 


E = E\ + E 2 + £ 3 , 
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( 40 ) 









where E\, E 2 are defined in equations f l2Jl) and 



f - Yj (w - w 7r dia S 




+ 2 Q,) 


(w - if/e) dx dz- 


(41) 


Then, E satisfies 


dE 

dt 


- f ~t 


+w T diag 


<A/ diag 


fld (0p + Q|) 


fixity 


+ 2 a) 


*Pt 


dp 




+ 2 £ 2 ;) 




► dxdz . 


(42) 


The proof of the proposition's similar to the proof of the proposition!]] and will not be repeated 
here. Proposition [3] calls the following comments: 

• only £3 and the time evolution of E are modified by the diffusive approximation; 

• the abscissae 0 l £ are always positive, as explained in § 13.41 but not necessarily the weights 
a l £ . Consequently, in the general case, we cannot say that the Biot-DA model is stable. 
However, in the particular case where the coefficients 6 l v a l £ are all positive, E is an energy, 
and yy < 0: the Biot-DA model is therefore stable in this case. 

Proposition 4. Let us assume that the abscissae 0 l £ have been sorted in increasing order 

0* <$,<•■■ <4, f = 1,3, (43) 

and that the coefficients 0 l £ , a l £ of the diffusive approximation 02l) are positive. Then zero is an 
eigenvalue with multiplicity 6 of S. Moreover, the 2 N + 2 non-zero eigenvalues ofS (denoted s l £ , 
€ = 1 , • • •, A + 1 ) are real positive, and satisfy 

0 < s\ < e\ + Q f < ■ ■ ■ < S 1 N < e l N + Eli < s l N+v i = 1,3. (44) 


Proposition [4] is proven in Appendix C As we will see in § |H the proposition [4] ensures the 
stability of the numerical method. Positivity of quadrature abscissae and weights is again the 
fundamental hypothesis. 


3.4. Determining the Biot-DA parameters 

For the sake of clarity, the space coordinates and the subscripts due to the anisotropy are 
omitted. The quadrature coefficients aim to approximate improper integrals of the form 


(d + n) l, 2 w(t) 


1 f°° 1 ^ 

- I — ff) do ^ V a ( !p{t. Of), 
n Jo y/0 
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(45) 
















Moreover, the positivity of the quadrature coefficients is crucial for the stability of the Biot-DA 
model and its numerical implementation, as shown in propositions [3] and |4] Two approaches can 
be employed for this purpose. While the most usual one is based on orthogonal polynomials, the 
second approach is associated with an optimization procedure applied to the viscous operators 

dm 


3.4.1. Gaussian quadratures 

Various orthogonal polynomials exist to evaluate the improper integral (l45l) . The first method, 
proposed in ill , is to use the Gauss-Laguerre quadrature formula, which approximates improper 
integrals over M + . Slow convergence of this method is explained and corrected in ES.it consists 
in replacing the Gauss-Laguerre quadrature by a Gauss-Jacobi quadrature, more suitable for 
functions which decrease algebraically. A last improvement, proposed in d, consists in using 
a modified Gauss-Jacobi quadrature formula, recasting the improper integral (l45l ) as 


1 

7 T 


r°° l 

I — ifr(6) d6 = 

Jo yfe 


1 

71 


X +l i ts 

(1 - ~GY( 1 + ~fffm 
1 71 €=\ 


(46) 


with the modified memory variable ij/ defined as 


m = 


_4_ 

(i -oy-\i + oy +3 


i + o 

1-6 





(47) 


The abscissae Of, which are the zeros of the Gauss-Jacobi polynomials, and the weights a f can 
be computed by standard routines [3b|]. In d, the author proves that for fractional derivatives 
of order 1 /2, the optimal coefficients to use are y = 1 and 6=1. The coefficients of the diffusive 
approximation Of and af © are therefore related to the coefficients Of and df (l46l) by 


61 = 



1 4 at 

7i (1 -§,)(!+ g,) 3 ' 


By construction, they are strictly positive. 


(48) 


3.4.2. Optimization procedures 

In 013, we proposed a different method to determine the coefficients Ot and at of the 
diffusive approximation (l45l) . This method is based on the frequency expressions of the viscous 
operators and takes into account the frequency content of the source. Our requirement is therefore 
to approximate the viscous operator F jkd (oj ) by F da (oj) (l38t in the frequency range of interest 
1 = [^min, ^>max] ? centered on the central angular frequency of the source. This leads to the 
minimization of the quantity^ 2 with respect to the abcissae Of and to the weights at 


X 


2 


z 

k= 1 


F DA (0J k ) 

~ 1 

2 K 

= v 

F JKD (oj k ) 


/ j 
k= 1 


N 

€=1 


(£l + j(o k ) 1/2 

Of + O + j (Ok 


2 


- 1 


(49) 


where the angular frequencies Of are distributed linearly in / on a logarithmic scale of K points 


~ ^min 


( ^max 
^min 


k -1 
K -1 


k = \ -K. 
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In I 26 L l27ll . the abcissae 0 £ were arbitrarily put linearly on a logarithmic scale, as (l5Qb . Only 
the weights a 1 were optimized with a linear least-squares minimization procedure of d49b . Some 
negative weights were obtained, which represents a major drawback, at least theoretically, since 
the stability of the Biot-DA model can not be guaranteed. 

To remove this drawback and improve the minimization of^ 2 , a nonlinear constrained opti¬ 
mization is developed, where both the abcissae and the weights are optimized. The coefficients Of 
and at are now constrained to be positive. An additional constraint 6^ < 6 max is also introduced 
to ensure the computational accuracy in the forthcoming numerical method (§ |4]). Setting 

9 ( = (e' e ) 2 , a e = (a' e ) 2 , (51) 


the number of constraints decreases from 3 N to N leading to the following minimization prob¬ 
lem: 


minx 2 , O’e ^ V^max- (52) 

((V n’ \ t T 


The constrained minimization problem (l52l) is nonlinear and non-quadratic with respect to ab¬ 
scissae Q' t To solve it, we implement the program SolvOpt l37[ 3ji], used in viscoelasticity & 
Since this Shor’s algorithm is iterative, it requires an initial estimate 0'°, a '° of the coefficients 
which satisfies the constraints of the minimization problem (l52b . For this purpose, 6® and a® 
are initialized with the method based on the modified Gauss-Jacobi quadrature formula (l48l ) . 
Different initial guess have been used, derived from Gaus-Legendre and Gauss-Jacobi methods, 
leading to the same final coefficients 6i and ct£. 

In what follows, we always use the parameters <jo m j n = u>o/10, c o max = 106l)o, 0 m ax = 

100 (j) o, K - 2 A, where coq = 2 n /o is the central angular frequency of the source. 


3.4.3. Discussion 

To compare the quadrature methods presented in § !3.4.1l and l3.4.2l we first define the error 
of model s mo d as 


^mod 


F da (oj) l 

| Z*<^max 

F da {cj) 1 

2 

daj 

F jkd (oj) 

L 2 \^ w min 

F JKD (oS) 

/ 


(53) 


The variation of s mo d in terms of the number A of memory variables, for /o = 200 kHz and 
f c = 3.84 kHz, is represented on figure [3]-a. The Gauss-Jacobi method converges very slowly, 
and the error is always larger than 1 % even for N = 50. Moreover, for values of N ^ 10, 
the error is always larger than 60 %. For both the linear and the nonlinear optimizations, the 
errors decrease rapidly with N. Nevertheless, the nonlinear procedure outperforms the results 
obtained in the linear case. For N = 8 for instance, the relative error of the nonlinear optimization 
(£mod - 7.16 10 -3 %) is 514 times smaller than the error of the linear optimization (s mo d - 3.68 
%). For larger values of A, the system is poorly conditioned and the order of convergence 
deteriorates; in practice, this is not penalizing since large values of N are not used. An example 
of a priori parametric determination of N in terms of both the frequency range and the desired 
accuracy is also given in figure [3]-b for the nonlinear procedure. The case N = 0 corresponds to 
the Biot-LF model. 

It is also important to compare the influence of the quadrature coefficient on the physical 
observables. For that purpose, we represent on figure 0 the phase velocity and the attenuation 
of the slow wave of the Biot-DA model, obtained with the different quadrature methods. As 
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(a) (b) 



Figure 3: (a): relative error s mo d in terms of N for both the modified Gauss-Jacobi quadrature and the nonlinear con¬ 
strained optimization, (b): required values of N in terms of fo/f c i and the required accuracy s moc i for the nonlinear 
optimization. 


expected, the results given by the Gauss-Jacobi method are extremely poor. On the contrary, 
the linear and non-linear procedures are able to represent very accurately the variations of these 
quantities on the considered range of frequencies, even for the small values N = 3. Based on 
these results and the positivity requirement, the nonlinear constrained optimization is therefore 
considered as the better way to determine the coefficients of the diffusive approximation. This 
method is used in all what follows. 

(a) (b) 




Figure 4: phase velocity (a) and attenuation (b) of the slow quasi-compressional wave. Comparison between the Biot-DA 
model and the Biot-JKD model for N = 3. 


4. Numerical modeling 

4.1. Splitting 

In order to integrate the Biot-DA system (l37l) . a uniform grid is introduced, with mesh size 
A x, f Az and time step A t. The approximation of the exact solution Ufa = i Ax, zj — j Az, t n = 
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nAt) is denoted by with 0 < i < N x , 0 < j < N z . If Ax = Az, a straightforward 

discretization of (1771) by an explicit time scheme typically leads to the following condition of 
stability 


( 


At ^ min 


Ax 


max c , 

(pe[0,n/2] PI 


(<fiY R(S) 


(54) 


where R(S ) is the spectral radius of S', and T > 0 is obtained by a Von-Neumann analysis when 
S = 0. The first term of (l54l h which depends of the propagation matrices A and B , is the 
classical CFL condition. The second term of (l54t depends only on the diffusive matrix S. From 
proposition @1 we deduce that the spectral radius of S satisfies 


R(S ) > max (Ol + kl\,6] + D 3 ) 


(55) 


if the coefficients 6 l { and ^ of the diffusive approximation are positive. With highly dissipative 
fluids, the second term of can be so small that numerical computations are intractable. 

A more efficient strategy is adopted here, based on the second-order Strang splitting Si. it 
consists in splitting the original system (f3Tb into a propagative part 


au A du 
+ ^ 

ot ox 


dU ^ 
+ B— = 0, 

oz 


(H p ) 


(56) 


and a diffusive part with forcing 


a u 

dt 


-SU + G , 


(H d ) 


(57) 


where H p and H d are the operators associated with each part. One solves alternatively the 
propagative part and the diffusive part: 


At 


U n+l = H d \t n+U — o H p (At) o H d lt n , — \ U n 


At 


(58) 


The discrete operator H p associated with the propagative part ([56b is an ADER 4 (Arbitrary 
DERivatives) scheme S3i. This scheme is fourth-order accurate in space and time, is dispersive 
of order 4 and dissipative of order 6 H42I1 . and has a stability limit T = 1. On Cartesian grids, 
ADER 4 amounts to a fourth-order Lax-Wendroff scheme. A general expression of the ADER 
scheme, together with its numerical analysis, can be found in the section 4-3 of the thesis 0. 

The solution of (l57l) is given by 




f 

J to 


'?o+A?/2 


-S (to+At/2-T) 


G(t) dr. 


(59) 


^e~ s ^U(t 0 )-(I-e- s T) S-'Gitk), 


with k = n or n + 1. The exponential matrix e~ SAt/2 is computed numerically using the (6,6) 
Pade approximation in the ’’scaling and squaring method” 114411 . Proposition [4] ensures that the 
numerical integration of the diffusive step ([57b is unconditionally stable 0]. Without forcing, 
i.e. G = 0, the integration of the diffusive part (l57b is exact. 
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The full algorithm is therefore stable under the optimum CFL condition of stability 


At = T 


Ax 


max 

^e[0,7r/2] 


pf 


(</>)’ 


1 , 


( 60 ) 


which is always independent of the Biot-DA model coefficients. Since the matrices A, B and S 
do not commute, the order of convergence decreases from 4 to 2. Using a fourth-order ADER 
scheme is nevertheless advantageous, compared with the second-order Lax-Wendroff scheme: 
the stability limit is improved, and numerical artifacts (dispersion, attenuation, anisotropy) are 
greatly reduced. 


4.2. Immersed interface method 

Let us consider two transversely isotropic homogeneous poroelastic media Do and Di sep¬ 
arated by a stationary interface T, as shown in figure [TJ The governing equations 071) in each 
medium have to be completed by a set of jump conditions. The simple case of perfect bonding 
and perfect hydraulic contact along T is considered here, modeled by the jump conditions JH: 

[v s ] = 0, [ w.n] = 0, [cr.w] = [ p] = 0. (61) 

The discretization of the interface conditions requires special care. A straightforward stair-step 
representation of interfaces introduces first-order geometrical errors and yields spurious numeri¬ 
cal diffractions. In addition, the jump conditions (l6Tb are not enforced numerically if no special 
treatment is applied. Lastly, the smoothness requirements to solve (l56b are not satisfied, decreas¬ 
ing the convergence rate of the ADER scheme. 

To remove these drawbacks while maintaining the efficiency of Cartesian grid methods, im¬ 
mersed interface methods constitute a possible strategy |4(| [47l [241] . The latter studies can be 
consulted for a detailed description of this method. The basic principle is as follows: at the ir¬ 
regular nodes where the ADER scheme crosses an interface, modified values of the solution are 
used on the other side of the interface instead of the usual numerical values. 

Calculating these modified values is a complex task involving high-order derivation of jump 
conditions <ED, high-order derivation of the Beltrami-Michell equation ^ and algebraic ma¬ 
nipulation, such as singular value decompositions. All these time consuming procedures can 
be carried out during a preprocessing stage and only small matrix-vector multiplications need 
to be performed during the simulation. After optimizing the code, the extra CPU cost can be 
practically negligible, i.e. lower than 1% of that required by the time-marching procedure. 

Compared with § 3-3 of 0, the modifications induced by anisotropy concern 

• step 1: the derivation of the jump conditions, 

• step 2: the derivation of the Beltrami-Michell equation. 

These modifications are tedious and hence will not be repeated here. They are deduced from the 
new expressions ([8]) and (l34l ) . 

4.3. Summary of the algorithm 

The numerical method can be summed up as follows: 

1. pre-processing step 
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• diffusive coefficients: initialisation (l48k nonlinear optimisation (l5Tl) - (l52t 

• numerical scheme: ADER matrices for (l56k exponential of the diffusive matrix (l59l) 

• immersed interface method: detection of irregular points, computation of extrapola¬ 
tion matrices 

2 . time iterations 

• immersed interface method: computation of modified values near interfaces 

• diffusive half-step Hj (l59l) 

• propagative step H p (l56l) , using modified values near interfaces 

• diffusive half-step H c j (l59t 


5. Numerical experiments 

Configuration 

In order to demonstrate the ability of the present method to be applied to a wide range of 
applications, the numerical tests will be run on two different transversely isotropic porous media. 
The medium Do is composed of thin layers of epoxy and glass, strongly anisotropic if the wave¬ 
lengths are large compared to the thickness of the layers [20]. The medium Di is water saturated 
Berea sandstone, which is sedimentary rock commonly encountered in petroleum eng ineering. 
The grains are predominantly sand sized and composed of quartz bonded by silica 120.14811. 

The values of the physical parameters are given in table [T] The viscous characteristic lengths 
Ai and A 3 are obtained by setting the Pride numbers P\ = P 3 = 0.5. We also report in these 
tables some useful values, such as phase velocities, critical frequencies, and quadrature param¬ 
eters computed for each media. The central frequency of the source is /o = 200 kHz, and the 
quadrature coefficients 6 l v a l v i = 1, 3, are determined by nonlinear constrained optimization with 
N = 3 memory variables. The error of model s mo d (l53l) is also given. We note that the transition 
frequencies f c \ and / c3 are the same for both Do and Di. In this particular case, the coefficients of 
the diffusive approximation are therefore also the same. In all the numerical simulations, the time 
step is computed from the physical parameters of the media through relations (l60l) , setting the 
CFL number T = 0.95. The numerical experiments are performed on an Intel Core i7 processor 
at 2.80 GHz. 

In the first test and the third test, the computational domain [-0.15,0.15] 2 m is discretized 
with N x - N z = 2250 grid nodes in each direction, which amounts to 20 points per slow com- 
pressional wavelength in Do. In the other tests, the computational domain [-0.1,0.1 ] 2 m is 
discretized with N x = N z = 1500, which amounts also to 20 points per slow compressional 
wavelength in Do and in Di. 

Test 1: homogeneous medium 

In the first test, the homogeneous medium Do (table [I]) is excited by a source point located at 
(0 m, 0 m). The only non-null component of the forcing F (l36l) is Go- n = git) h(x , z), where g(t) 
is a Ricker signal of central frequency /o and of time-shift to = 1 //o = 10 -5 s: 

I (2 7T 2 /o (t - to) 2 - 1 ) exp (-n 2 / 0 2 (t-to) 2 ) if 0 < t < 2 1 0 , 

Hit) = | (62) 

[ 0 otherwise, 
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Parameters 

Do 

Oi 

Saturating fluid 

p f (kg/m 3 ) 

1040 

1040 


rj (Pa.s) 

io - 3 

10“ 3 


K f (GPa) 

2.5 

2.5 

Grain 

p s (kg/m 3 ) 

1815 

2500 


K s (GPa) 

40 

80 

Matrix 

0 

0.2 

0.2 


r, 

2 

2 


r 3 

3.6 

3.6 


ki (m 2 ) 

6 .10- 13 

6 . IO" 13 


k 3 (m 2 ) 

10“ 13 

10“ 13 


c \i (GPa) 

39.4 

71.8 


ci2 (GPa) 

1 

3.2 


cu (GPa) 

5.8 

1.2 


c 33 (GPa) 

13.1 

53.4 


c 55 (GPa) 

3 

26.1 


At (m) 

6.93 10“ 6 

2.19 IO" 7 


A 3 (m) 

3.79 10“ 6 

1.20 IO" 7 

Dispersion 

c p f (0 ) (m/s) 

5244.40 

6004.31 


Cpf(fo, 0) kHz (m/s) 

5227.10 

5988.50 


c“/*r/2) (m/s) 

3583.24 

5256.03 


Cpf{fo,7i/2 ) (m/s) 

3581.42 

5245.84 


(0) (m/s) 

975.02 

1026.45 


Cpsifo, 0 ) (m/s) 

901.15 

949.33 


cJ>/2) (m/s) 

604.41 

745.59 


Cp S (fo,n/2) (m/s) 

534.88 

661.32 


cj°(0) (m/s) 

1368.36 

3484.00 


c,(/ 0 ,0) (m/s) 

1361.22 

3470.45 


c“(tt/2 ) (m/s) 

1388.53 

3522.07 


c 5 (/ 0 ,7t/ 2) (m/s) 

1381.07 

3508.05 


/cl (Hz) 

2.55 10 4 

2.55 10 4 


/c3 (Hz) 

8.50 10 4 

8.50 10 4 

Optimization 

0\ (rad/s) 

1.6410 5 

1.6410 5 


6\ (rad/s) 

2.80 10 6 

2.80 10 6 


6 \ (rad/s) 

3.58 10 7 

3.58 10 7 


aj (rad 1/2 /s 1/2 ) 

5.58 10 2 

5.58 10 2 


(rad 1/2 /s 1/2 ) 

1.21 10 3 

1.21 10 3 


(rad 1/2 /s 1/2 ) 

7.32 10 3 

7.32 10 3 



1.61 

1.61 


0 3 (rad/s) 

3.14 10 s 

3.14 10 5 


6» 3 (rad/s) 

5.06 10 7 

5.06 10 7 


6> 3 (rad/s) 

4.50 10 6 

4.50 10 6 


a 3 (rad 1/2 /s 1/2 ) 

7.57 10 2 

7.57 10 2 


a 3 (rad 1/2 /s 1/2 ) 

8.79 10 3 

8.79 10 3 


a 3 (rad 1/2 /s 1/2 ) 

1.38 10 3 

1.38 10 3 



0.53 

0.53 


Table 1: Physical parameters of the transversely isotropic media used in the numerical experiments. The phase velocities 
Cpf(fo,(f), Cp S (fo, <p) and c s (/o, p) are computed at / = /o = 200 kHz when the wavevector k makes an angle ip with the 
horizontal x-axis, and c^(cp), c™ s {<p), c™{<p) denote the high-frequency limit of the phases velocities. 
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and h(x,z) is a truncated Gaussian centered at point (0,0), of radius Ro = 6.56 10 3 m and 
E = 3.28 10“ 3 m: 


h(x, z) = 



X2+Z 2 \ 

E 2 ) 


if 0 < x 2 + z 2 < R 2 0 , 


0 


otherwise. 


(63) 


The time step is At = 2.41 10 -8 s. We use a truncated Gaussian for h(x,z) rather than a Dirac 
distribution to avoid spurious numerical artifacts localized around the source point. This source 
generates cylindrical waves of all types: fast and slow quasi-compressional waves and quasi¬ 
shear waves, which are denoted by qPf , qP s and qS , respectively, in figure 0 The three waves 
are observed in the pressure field. Contrary to the isotropic case, where the pressure of the shear 
wave is null, pressure is visible in the qS wave. 

A comparison is proposed with the theoretical wavefronts deduced from the dispersion anal¬ 
ysis (section [276b and the resolution of d29b . They are denoted by a black dotted line in figure [5] 
It is observed that the computed waves are well positionned at the final instant t\ ^ 2.72 10 -5 s 
(corresponding to 1125 time steps). No special care is applied to simulate outgoing waves (with 
PML, for instance), since the simulation is stopped before the waves have reached the edges of 
the computational domain. The cusp of the shear wave is seen in the numerical solution. 


P 


zoom 



x (m) 


Figure 5: test 1. Fast and slow quasi-compressional waves, respectively qPf and qP s , and quasi-shear wave qS emitted 
by a source point at (0 m, 0 m). Pressure at t\ 2.72 10 -5 s. 

Test 2: diffraction of a plane wave by a plane interface 

In the second test, the source is a plane fast compressional wave traveling in the positive 
direction of the v-axis, whose wavevector k is parallel to the direction of propagation. Its time 
evolution is the same Ricker signal as in the first test ([62b . We use periodic boundary conditions 
at the top and at the bottom of the domain. The validity of the method is checked in the particu¬ 
lar case of heterogeneous transversely isotropic media, where a semi-analytical solution can be 
obtained easily. The media Qo and Qi are separated by a vertical wave plane interface at v = 0 
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m. The incident Pf -wave (Jpf) propagates in the medium Qi. The time step is At = 2.11 10 -8 
s. The figure [6] shows a snapshot of the pressure at t\ - 1.48 10 -5 s (corresponding to 750 time 
steps), on the whole computational domain. The reflected fast and slow quasi-compressional 
waves, denoted respectively Rpf and Rp s , propagate in the medium Di; and the transmitted 
fast and slow quasi-compressional waves, denoted respectively Tpf and Tp s , propagate in the 
medium Do. In this case, we compute the exact solution of Biot-DA thanks to Fourier tools and 


(a) 


(b) 




Figure 6: test 2. Snapshot of pressure at initial time (a) and at t\ =* 1.48 10 5 s (b). The plane interface is denoted by a 
straight black line, separating T>] (on the left) and Qq (on the right). 


poroelastic equations; a general overview of the analytical solution is given in the Appendix D 


The figure [7] shows the excellent agreement between the analytical and the numerical values of 
the pressure along the line z = 0 m. Despite the relative simplicity of this configuration (ID 
evolution of the waves and lack of shear waves), it can be viewed as a validation of the numerical 
method which is fully 2D whatever the geometrical setting. 


Test 3: source point and plane interface or sinusoidal interface 

In the previous test, the configuration was fully ID, but more complex geometries can be han¬ 
dled on a Cartesian grid thanks to the immersed interface method. As an example, the media Do 
and Di are separated by a plane interface with slope 15 degree with the horizontal v-axis, passing 
through the point (0 m, -0.004 m). The homogeneous medium Di is excited by the source point 
described in test 1. This source emits cylindrical waves which interact with the medium Do- The 
time step is At = 2.11 10 -8 s. Snapshot of the pressure at time t - 2.53 10 -5 s (corresponding 
to 1200 time steps) and the time evolution of the pressure at receivers Rq (0.042 m, 0.0068 m) in 
Do and R\ (0.048 m, 0.0071 m) in Di are represented on figured 

The plane interface can be easily replaced, for instance by a sinusoidal interface of equation 

- sin 6(x- x s ) + cos 6(z- z s ) ~ A s sin (oj s (cos 6(x- x s ) + sin 6{z~ z s ))) - 0, (64) 

with x s = 0 m, z s = -0.027 m, A s = 0.01 m, oj s = 50 n rad/s, 6 = tt/ 12 rad. Snapshot of the 
pressure and the time evolution of the pressure at receivers Rq (0.036 m, -0.031 m) in Do and R\ 
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Pressure 


Zoom on the slow compressional waves 




Figure 7: test 2. Pressure along the line z = 0 m; vertical line denotes the interface. Comparison between the numerical 
values (circle) and the analytical values (solid line) of p at t\ ^ 1.48 10 -5 s. 
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Figure 8: test 3. Snapshot of pressure and at t\ ^ 2.53 10 -5 s (a). The plane interface separating the media F>o and £l\ 
is denoted by a straight black line. Time evolution of the pressure (b) at receiver Ro in Qq (red) and at receiver R\ in 
(blue). 
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(0.036 m, -0.014 m) in Di are represented on figured In both cases, no spurious diffraction is 



Figure 9: test 3. Snapshot of pressure and at t\ - 2.53 10 5 s (a). The sinusoidal interface separating the media Qq and 
i is denoted by a straight black line. Time evolution of the pressure (b) at receiver Rq in F>o (red) and at receiver R\ in 
Qi (blue). 

induced by the stair-step representation of the interface, thanks to the immersed interface method. 
Moreover, classical waves conversions and scattering phenomena are observed: reflected, trans¬ 
mitted and Stoneley waves. The shape of the transmitted waves, not circular, illustrates the strong 
anisotropy of the medium Do. 

Test 4: multiple ellipsoidal scatterers 

To illustrate the ability of the proposed numerical strategy to handle complex geometries, 
200 ellipsoidal scatterers of medium Di, with major and minor radii of 0.025 m and 0.02 m, are 
randomly distributed in a matrix of medium Do- The computational domain is [-0.8,0.8] 2 m, 
hence the surfacic concentration of scatterers is 25 %. A uniform distribution of scatterers is 
used. The source is the same plane fast compressional wave than is the second test, and we use 
periodic boundary conditions at the top and at the bottom of the domain. 

The time step is At = 3.37 1CT 7 s. The pressure is represented at the initial time and at 
time t\ - 1.43 10~ 4 s (corresponding to 425 time steps) on figure This simulation has taken 
approximately 11.5 h of preprocessing and 8.5 h of time-stepping. Similar numerical experiments 
are also performed for a surfacic concentration of scatterers of 10 % and 15 %. 

At each time step, the components of £/"• are stored inside the subdomain containing the 
inclusions. For this purpose, a uniform network consisting of A/ = 800 lines and N c = 25 
columns of receivers is put in the subdomain. The position of the receivers is given by ( x^zj ), 
where i = 0, • • ■, N c - 1 and j = 0, • • •, A/ - 1. The field £/"• recorded on each array (each line 
of receivers), represented on figureQjJa, corresponds to a field propagating along one horizontal 
line of receivers. A main wave train is clearly visible, followed by a coda. Summing the time 
histories of these N c arrays gives a coherent field propagating in the v direction: 
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Figure 10: test 4. Multiple scattering in random media. Snapshot of p at the initial time (a) and at time t\ « 1.43 10 4 s. 
The matrix is whereas the 200 scatterers are F>i. 


On the coherent seismogram thus obtained, represented on figureQjJb, the coda has disappeared, 
and the main wave train behaves like a plane wave propagating in a homogeneous (but disper¬ 
sive and attenuating) medium. The coherent phase velocity c(oj), represented in figure [12]- a, is 
computed by applying a p - oj transform to the space-time data on the coherent field (l65l ) . where 
p = l/c is the slowness of the waves [49, 50]. The horizontal lines represent a simple average of 
the phase velocities weighted by the concentration. The coherent attenuation a(co) is estimated 
from the decrease in the amplitude spectrum of the coherent field during the propagation of the 
waves, see[l2]-b. An error estimate is also deduced, represented in figure [12] by vertical lines. 


6. Conclusion 

An explicit finite-difference method has been developed here to simulate transient poroelastic 
waves in the full range of validity of the Biot-JKD model, which involves order 1 /2 fractional 
derivatives. A diffusive representation transforms the fractional derivatives, non-local in time, 
into a continuum of local problems, approximated by quadrature formulae. The Biot-JKD model 
is then replaced by an approximate Biot-DA model, much more tractable numerically. The co¬ 
efficients of the diffusive approximation are determined by a nonlinear constrained optimization 
procedure, leading to a small number of memory variables. The hyperbolic Biot-DA system of 
partial differential equations is discretized using various tools of scientific computing: Strang 
splitting, fourth-order ADER scheme, immersed interface method. It enables to treat efficiently 
and accurately the propagation of transient waves in transversely isotropic porous media. 

Some future lines of research are suggested: 

• Multiple scattering. Many theoretical methods of multiple scattering have been developed 
to determine the effective wavenumber of media with random scatterers; see for instance 
the Independent Scattering Approximation and the Waterman-Truell method & The 
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Figure 11: test 4. Incident plane qP s -wave in a medium with 25% inclusion concentration, (a): pressure recorded along 
an array, (b): coherent pressure obtained afer summation. 
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Figure 12: test 4. Effective phase velocity (a) and effective attenuation (b) at various inclusion concentrations. The 
vertical lines represents the error bars. The horizontal lines in (a) give the average phase velocity weighted by the 
concentration. 
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main drawback of these methods is that their validity is restricted to small concentrations 
of scatterers, typically less than 10 %. On the contrary, numerical methods do not suf¬ 
fer from such a limitation if suitable efforts are done. In particular, the errors due to the 
discretization (numerical dispersion, numerical dissipation, spurious diffractions on inter¬ 
faces, ...) must be much smaller than the physical quantities of interest. In l52ll . numerical 
simulations were used in the elastic case to estimate the accuracy of standard theoretical 
models, and also to show the improvement induced by recent models of multiple scattering 
153]. As shown in test 4 of §[5] the numerical tools presented here make possible a similar 
study poroelastic random media and comparisons with theoretical models lHH \5^. 

However, realistic configurations would involve approximately 1500 scatterers, and sizing 
of the experiments leads to N x xN z = 10000 2 , and 10000 time iterations are required. Con¬ 
sequently, the numerical method has to be parallelized, for instance by Message Passing 
Interface (MPI). 

• Thermic boundary-layer. In cases where the saturating fluid is a gas, the effects of thermal 
expansion of both pore fluid and the matrix have to be taken into account. In the HF 
regime, the thermal exchanges between fluid and solid phase occur in a small layer close 
to the surface of the pores. In this case, the dynamic thermal permeability is introduced 
[56], leading in the time-domain to an additional shifted fractional derivative of order 
1/2. The numerical method developed in this paper can be applied without difficulty by 
introducing additional memory variables. 

• Fractional derivatives in space. The Biot theory is very efficient to predict the macro¬ 
scopic behavior of long-wavelength sound propagation in porous medium with relatively 
simple microgeometries. However, it remains far to describe correctly the coarse-grained 
dynamics of the medium when the microgeometry of the porous medium become more 
complex, for instance fractal. For rigid-framed porous media permeated by a viscothermal 
fluid, a generalized macroscopic nonlocal theory of sound propagation has been developed 
to take into account not only temporal dispersion, but also spatial dispersion £ 7 ]. In this 
case, the coefficients depends on the frequency and on the wavenumber. In the space-time 
domain, it introduces not only time-fractional derivatives, but also space-fractional deriva¬ 
tives. Numerical modeling of space-fractional differential equations has been addressed 
by several authors |j58L 85 9 1. by using a Grunwald-Letnikov approximation. The diffusive 
approximation of such derivatives constitutes an interesting challenge. 
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Appendix A. Proof of proposition Q] 

The equation (llOab is multiplied by v s T and integrated 



(A.l) 
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The first term in (1A.11) is written 

r T d v s , , d 1 P Til 

pv s -—dxdz=-r- pv s V s dxdz. 

Jr 2 dt dt 2 J r2 

Integrating by part and using ©, we obtain 

- f v s T (V.o-)dxdz = f cr T dxdz = f cr T (c~ x ^ - C~ l p ^ 

Jr 2 Jr 2 at J R 2 \ dt dt 

= —r — f cr T C _1 crdxdz + f cr T C _1 p'J- dxdz, 
dt 2 Jr2 J r2 dt 

--r\ f (cr T C~ l cr + 2 a T C~ l P p) dx dz - f (-r- ) C _1 pp dxdz. 
dt 2 Jr2 ' > J R 2\dtJ 

Equation (llObb is multiplied by w T and integrated 

X f r r d V g T i' / \ ^ ^ T T-r 

\Pf w TT + w diag (p wi ) —- + w Vp 
2 { J dt dt 

+w T diag(^ ^ (Z) + Q*) 1/2 ) w} dxdz = 0. 

The second term in (1A.41) can be written 

r T dw dlC T 

w dia g(p wi ) —dxdz = — - w dia g(p wi )wdxdz. 

Jr 2 dt dt 2 J R2 

Integrating by part the third term of (1A.4K we obtain 
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Equations (fl9l) and dA. lb - dA.71 ) yield 


(A.2) 


(A.3) 


(A.4) 


(A.5) 


(A.6) 


1 f 1 , C T , dcr d 1 r T , 

= —- I —p 2 dxdz+ I P T C 1 —— pdxdz+ — - I p T C x pp 2 dxdz. 
dt 2 Jr2 ot J K 2 dt dt 2 J r2 

We add dA.lb and the three first terms of dA.4b . Using the symmetry of C, there remains 


(A.7) 


— (E| + £ 2 ) = - f f --w^diagf-—) tj/dOdxdz. (A.8) 

dt Jr 2 Jo 7rV0 \KiVtoi) 
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To calculate the right-hand side of (1A.81) . equation (121al) is multiplied by w T or i// T 

( T dllf T dw T T 

w T ~x~ -w T — +w T diag (6 + A) ij/-w T diag (A) w = 0 , 
ot ot 

V lu ~ 'P T7 + 'P diag + Qi) ^ - 'P diag ^ w = 0 

Ot ot 


Equation (IA.9I) can be written as 


d 1 


y diag (0 + 2 A) w = (w - i/f) 1 (w - i/f) 

ot 2 

+^r r diag (0 + A) + hA diag (A) w. 
Substituting (1A.101) in (IA.8I) leads to the relation (f24l) 
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(A.9) 


(A. 10) 


(A.l 1) 


’ W; VTT(0+ 2Qf)J 

It remains to prove that £ d22l) is a positive definite quadratic form. Concerning £i,we write 


pv/v s + *r r diag(p wi ) w + 2p/v/w = X T 7 # 1 X 1 + X 3 r # 3 X 3 , 


(A. 12) 


where 


Xj = (v,„- w,) 7 ’, 


//, = 


P Pf 

Pf Pwi , 


i= 1,3. 


(A.13) 


Taking <S; and Pi to denote the sum and the product of the eigenvalues of matrix Hi, we obtain 


( !P, = det Hi = pp wi -pj= Xi > 0, 
( Si = lx Hi = p+p w > 0 . 


(A.14) 


The eigenvalues of Hi are therefore positive. This proves that E\ is a positive definite quadratic 
form. The terms £’ 2 , £3 and are obviously positive definite quadratic form because the 
involved matrices are definite positive. □ 


Appendix B. Matrices of propagation and dissipation 


The matrices in (l37l) are 
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Appendix C. Proof of proposition 0| 

We denote P<g the change-of-basis matrix satisfying 

U = P S (U l9 U 3 ,cr, p) T , 

with 


(C.l) 

(C.2) 


Vi = (v si i // l N ) T , i = 1, 3. 

The matrix P<g is thus invertible, and the matrices S ( [Appendix B| ) and S$ = S'Pg are similar. 
The matrix S% writes 
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with (i = 1,3) 
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The characteristic polynomial of S is 

P S (s) = S 4 P S (s)P S (s), 


(C.5) 
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where Pg (s) denotes the characteristic polynomial of the matrix Si, i.e. Si(s ) = det(S; - sI^+2) 
with I N+ 2 the (N + 2)-identity matrix. This (N + 2)-determinant is expanded along the first 
column. The line I and the column J of the (N + l)-determinant thus obtained are denoted L/ 
and Cj, respectively (0 ^ I, J ^ N). The following algebraic manipulations are then performed 
successively: 


(0 U^U-U, €=1 ,---,N, 

(ii) C 0 <- C 0 n(^ + Q ~ s), 
e=\ 

(Hi) Co^Co-(s-£h) + €=1 

k=i 

k-H 

One deduces 

N N N 

P §i (s) = ~ x Qi(s) = S 2 + Q - s) + 7i s(s- £li) 2 a\ ]~[(4 + Q, - s). 

€=1 £=1 k= 1 

kU 


(C.6) 
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From equation (1C.7F one has (0) ^ 0 while (2/(0) + 0, therefore 0 is an eigenvalue of the 

matrix Si with multiplicity 1. In what follows, the positivity of the coefficients 6 l p a l £ of the 
diffusive approximation is used. In the limit s —> 0 + , then asymptotically 


P §i (s) ~■ -7,0, s 24 fl ( ^ + a) => sgn(f$(0 + )) = -1- (C.8) 

+ €=1 k= 1 \ * / 

kU 

Moreover, using (l43lh then at the quadrature abscissae one has for all € = 1 , • • •, TV 

N ( v 

= r, 4 (4 + n,)a‘ F[(4 - 4) => sgn(p^ (4 + .Q,)j = (-1/+ 1 . (C.9) 

k—\ 

Finally, the following limit holds 

PcO) ~ i-lfs"* 2 ^sgn(p, (+00)) = (-If. (C.10) 

s—»+oo \ / 

We introduce the following intervals 

4 =]ff N + Qi, +oo[, 4=]^,4 1+ Od, forT’ = 1, • • •,TV — 1, /j =10,6* + Q f ]. (C.ll) 

The real-valued continuous function Pg changes of sign on each interval I l r Consequently, 
according to the intermediate value theorem, Pg has at least one zero in each interval. Since Pg 
has at the most N + 1 distinct zeros in ]0, +oo[, we deduce that 3 ! s' £ e I\IP$ (s\) = 0, l - 
1, • • • , TV + 1. Using equation (1C.5F the characteristic polynomial of S (1C. 71) is therefore 


P S (S) = s 6 
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T 5 - 4) (■* - 4)> 


e=\ 


(C.12) 


which concludes the proof. □ 
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Appendix D. Semi-analytical solution of Test 2 


In test 2, we consider the interaction of a plane wave with a plane interface at normal inci¬ 
dence. This ID case can be solved semi-analytically, by using Fourier analysis and poroelastic 
constitutive equations. Note that no shear wave is involved here. The general overview of the 
algorithm is as follows: 


• writing the potentials of fluid and elastic motions in terms of the potentials of fast and 
slow compressional waves. To do so, diagonalize the vectorial Helmholtz decomposition 
of Biot equations H]. The coefficient 

((1 - <f>)p s +p f P(T - 1)) co 2 - U f + 2p-mj3 2 ) k 2 + iaj(f) 

Y(co) = - ---^---- (D.l) 

Pf(T/3-(p)(jj 2 -ioj(/)j3 -F(oj) 

K 

is introduced by the change of basis, where k is the wavenumber and F(co) is the JKD 
frequency correction; 

• deduce each field (velocity, stress, pressure) from the expressions of potentials and from 
the poroelastic constitutive equations; 

• for a given Fourier mode, the reflected and transmitted waves can be written 


U(jc, co) 


+ 1 
0 

±0(1 - Y) 
0 


— [Af + 2/j + fim(/)(Y - 1)) 
0 

^(Af+/3mc/>(Y- 1)) 
——m (J3 + (p(Y - 1)) 


i(ut-kx) 


g(u). 


(D.2) 


where g(co) is the Fourier transform of the source 


The sign ± depends whether the 


wave is a right-going incident or transmitted wave (+) or a left-going reflected wave (-); 

• the four reflected and transmitted slow and fast waves (ID. 21) must be multiplied by a reflec¬ 
tion or transmitted coefficient. The four coefficients are computed by applying the jump 
conditions and by solving the resulting linear system; 


an inverse Fourier transform yields an approximate time-domain solution. 
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